home *** CD-ROM | disk | FTP | other *** search
/ TeX 1995 July / TeX CD-ROM July 1995 (Disc 1)(Walnut Creek)(1995).ISO / graphics / gnuplot / contrib / byrne / congp2d3.for < prev    next >
Text File  |  1992-03-25  |  25KB  |  508 lines
  1. c     C  O N G P 2 D 3   last revised Feb 6, 1992
  2. c     This program is made available courtesy of Rose Byrne
  3. c     Feb 10, 1992.  Rose is a doctoral candidate in Civil Engineering
  4. c     at UVa.
  5. c
  6. c     Congp2d3.for is a revision of congp2dr.for.  Congp2d3.for
  7. c     extends the algorithm described in Roache, "Computational
  8. c     Fluid Mechanics", to an irregular, simply connected,
  9. c     possibly non-convex region.  If the region is non-convex, the
  10. c     indentation should be placed in such a way that there is no
  11. c     node j which is neither physically adjacent to node j-1
  12. c     nor along the side whose nodes are specified in nbot.
  13. c     The algorithm in Roache originally appeared in Sundberg,
  14. c     "Two Computerized Methods For Plotting Functions of Two
  15. c     Independent Variables". This algorithm involves searching each
  16. c     element for sides which contain points on a given contour.
  17. c     The points within the element are then connected, and the pen
  18. c     is lifted.  Even complicated results can be contoured
  19. c     WITHOUT MANUALLY REORDERING THE DATA POINTS.  In the rare cases
  20. c     when a contour intersects a cell in four places, the program
  21. c     prepares a data file which connects those four dots in all 6
  22. c     possible ways.  The user should plot the file first and determine
  23. c     from the context of the plot how those four points should be
  24. c     connected.  Inserting a blank line in the data file will cause
  25. c     gnuplot to pick up the pen.  Congp2d3.for finds the points and
  26. c     writes a data file for gnuplot with blank lines in each place
  27. c     where the pen should be lifted.  Congp2d3.for also lifts the pen
  28. c     in each place where the region is non-convex.  In each place
  29. c     where the region is non-convex, the contours must stop at the
  30. c     boundary of the region and pick up on the other side.
  31. c     Sundberg's algorithm as extended here treats contours on
  32. c     non-convex regions correctly without the user's intervention.
  33. c     However, the section of the code which identifies points on the
  34. c     boundary cannot distinguish between a rectangular indentation
  35. c     in the region and a trapezoidal indentation.  The points which
  36. c     are on the boundary of the region but are also in the same
  37. c     column as another point on the boundary will not be found and must
  38. c     be entered manually.  For most regions, this will be only a couple
  39. c     of points whose coordinates will be obvious.  Process all data
  40. c     files for a given region and then fix the data file for the
  41. c     boundary.
  42. c
  43. c
  44. c     Congp2d3.for assumes that the points to be contoured are
  45. c     numbered as they usually are for finite element methods:
  46. c     from the bottom to the top of each column, with each column
  47. c     beginning where the last one leaves off.  Congp2d3.for keeps
  48. c     track of which nodes are in each element using an array of
  49. c     the points at the bottom of the region.  The following
  50. c     is an example of how the points are numbered:
  51. c
  52. c           5-----10               23----28
  53. c           |     |                 |     |    33
  54. c           |     |                 |     |     |
  55. c           4-----9----14----18----22----27----32----37
  56. c           |     |     |     |     |     |     |     |    41
  57. c           |     |     |     |     |     |     |     |     |
  58. c           3-----8----13----17----21----26----31----36----40
  59. c           |     |     |     |     |     |     |     |     |   44
  60. c           |     |     |     |     |     |     |     |     |    |
  61. c           2-----7----12----16----20----25----30----35----39---43
  62. c           |     |     |     |     |     |     |     |     |    | 46
  63. c           |     |     |     |     |     |     |     |     |    |  |
  64. c           1-----6----11----15----19----24----29----34----38---42-45
  65. c
  66. c           the nodes on the bottom are 1,6,11,15,19,24,29,34,38,42,45
  67. c
  68. c     This program processes a certain kind of data file which cannot
  69. c     be contoured by gnuplot3.0 (even after re-ordering the points) so
  70. c     that gnuplot can construct a contour plot by drawing the contour
  71. c     lines and domain boundary on a 2-D grid using the 2-D version
  72. c     of gnuplot.  This program processes a data file such as that
  73. c     produced by a finite element or finite difference program which
  74. c     contains values of the dependent variable at each point of a 2-D
  75. c     possibly irregular grid on a possibly irregular domain.  The
  76. c     elements are not necessarily rectangular.  If triangular elements
  77. c     are used in the finite element model which generates the data,
  78. c     the side of a triangular element which is on the boundary will not
  79. c     be searched for contour crossings, but the elements on the
  80. c     interior of the region will be combined to make quadrilateral
  81. c     elements.  Non-rectangular quadrilaterals work fine.
  82. c
  83. c     This program reads in the three arrays containing the three
  84. c     coordinates (dependent variable and two independent variables)
  85. c     of points representing the data which the user would
  86. c     like to have contoured.  These points may correspond to the points
  87. c     of any finite element or finite difference grid.  The arrays are
  88. c     assumed to be ordered in the order that the grid points are
  89. c     numbered. It is very helpful if the order in which the grid points
  90. c     are numbered is nowhere near parallel to the expected direction of
  91. c     the contour lines.  The program uses linear interpolation to
  92. c     construct level curves at equally spaced intervals:  the
  93. c     curves which will form the contours.  The coordinates of these
  94. c     curves are written to a file in the form which gnuplot requires.
  95. c     The curves will not contain the same number of points, so a
  96. c     contour plot cannot be constructed using the contour option of
  97. c     the 3-D version of gnuplot.  A contour plot can, however, be
  98. c     constructed using the 2-D version of gnuplot by drawing the
  99. c     independent variable coordinates of each contour line on the
  100. c     same 2-D grid.
  101. c
  102.       real zero
  103.       integer npts,nel,nlvl,lbot,mxlen
  104. c
  105. c     define parameters so dimension statements don't have to be
  106. c     changed each time the data is changed, yet arrays are no
  107. c     larger than they need to be.  The parameter statements must be
  108. c     changed each time the 2-D grid to be contoured or the number of
  109. c     contours desired is changed.
  110. c     npts = the number of points in the 2-D grid.
  111. c     nel = the number of elements in the 2-D grid.
  112. c     nlvl = the number of equally spaced level curves (contours)
  113. c            to be drawn.
  114. c     lbot = the number of points along the bottom (or lhs) of the grid.
  115. c     mxlen = the maximum number of points allowed in a contour.
  116.       parameter(npts=595)
  117.       parameter(nel=528)
  118.       parameter(nlvl=15)
  119.       parameter(lbot=52)
  120.       parameter(mxlen=400)
  121. c
  122. c     The program uses the following arrays:
  123. c     x = the array of the independent variable in one direction.
  124. c     z = the array of the independent variable in the other direction.
  125. c     h = the array of the dependent variable.
  126. c     xcont = the array of x coordinates of points on the contours.
  127. c            The contents, and the contents of zcont and the
  128. c            appropriate value of h for each contour will be
  129. c            written to the output file for processing by gnuplot.
  130. c     zcont = the array of z coordinates of points on the contours.
  131. c     len = the number of contour points in each contour line, not
  132. c           counting repetitions for Sundberg's connect the dots
  133. c           algorithm but counting each dot in each element in which
  134. c           it is found.
  135. c     value = the value of the dependent variable along a given contour.
  136. c     nbot = the array of points along the bottom of the grid (if points
  137. c            are numbered bottom to top).  If points are numbered left
  138. c            to right, this array contains the points on the left-hand
  139. c            edge. npts+1 is always added to the end of the array
  140. c     ncell = the number of contour points found in each element
  141. c     last = a running tally of the number of points in this contour
  142. c            level including repetitions for Sundber's algorithm.
  143.       real x(npts),z(npts),h(npts),xcont(nlvl,mxlen)
  144.       real value(nlvl), zcont(nlvl,mxlen)
  145.       integer ncell(nlvl,nel), last(nlvl), len(nlvl), nbot(0:lbot)
  146. c
  147. c     call subroutine readin to read from the screen the names of the
  148. c     input and output files and variables not included in the files,
  149. c     to open the necessary files, and to read the data from the input
  150. c     files.
  151.       call readin(zero,npts,nlvl,lbot,h,x,z,nbot,value,len)
  152. c
  153. c     call subroutine scan to scan all points in the array, looking for
  154. c     places where a contour value should be.
  155.       call scan(h,x,z,nbot,value,len,ncell,xcont,zcont,
  156.      $     npts,nel,nlvl,lbot,mxlen,zero)
  157. c
  158. c     call subroutine wrcon to write the information about the contours
  159. c     to a file in the correct format for gnuplot.
  160.       call wrcon(ncell,xcont,zcont,last,nlvl,nel,mxlen)
  161. c
  162. c     call subroutine wrbdy to write the information about the
  163. c     coordinates of the points on the boundary to another file
  164. c     for use by gnuplot.
  165.       call wrbdy(x,z,nbot,npts,lbot)
  166. c
  167. c     call subroutine wrlabl to write the labels for the contours
  168. c     to a third file for use by gnuplot.
  169.       call wrlabl(xcont,zcont,value,len,nlvl,mxlen)
  170. c
  171. c     To plot the data using gnuplot, first enter gnuplot by typing
  172. c     gnuplot with the dafault drive set to the drive that gnuplot
  173. c     is on.  Set any desired titles, axis labels, and output device
  174. c     (with appropriate options).  See the gnuplot documentation for
  175. c     details.  Then type the following commands:
  176. c         load "labels"
  177. c         plot "outcon", "outbnd"
  178. c     Insert the appropriate file names for labels, outcon, and outbnd
  179. c     above.  The quotation marks are part of the commands.
  180. c     See the gnuplot documentation for how to set axis ranges and
  181. c     many other features of the plot.
  182.       stop
  183.       end
  184. c
  185.       subroutine readin(zero,npts,nlvl,lbot,h,x,z,nbot,
  186.      $     value,len)
  187.       character*16 incont,outcon,rdnbot,outbnd,labels
  188.       integer npts, nlvl, lbot, iflag
  189.       real zero, to, from, by
  190.       real h(npts), x(npts), z(npts), value(nlvl)
  191.       integer nbot(0:lbot), len(nlvl)
  192. c     prompt user for input and output file names
  193.       write(*,*) 'Welcome to congp2d3, a preprocessor for gnuplot.'
  194.       write(*,*) 'This program will take the output from your 2-D'
  195.       write(*,*) 'finite element or finite difference program'
  196.       write(*,*) '(or a similar data file) and put it in a form'
  197.       write(*,*) 'in which it can be used by gnuplot.'
  198.       write(*,*)
  199.       write(*,*) 'Please type the name of the file with your results'
  200.       read(*,'(a16)') incont
  201.       write(*,*) 'Please type the name of the file with the array of'
  202.       write(*,*) 'points along the bottom of the grid (for points'
  203.       write(*,*) 'numbered bottom to top'
  204.       read(*,'(a16)') rdnbot
  205.       write(*,*) 'Thank you.  Now please type the name of the file to'
  206.       write(*,*) 'store the contour data for use by gnuplot.'
  207.       read(*,'(a16)') outcon
  208.       write(*,*) 'Thank you.  Now please type the name of the file to'
  209.       write(*,*) 'store the boundary data for use by gnuplot.'
  210.       read(*,'(a16)') outbnd
  211.       write(*,*) 'Thank you.  Now please type the name of the file to'
  212.       write(*,*) 'store the label definitions for gnuplot'
  213.       read(*,'(a16)') labels
  214.       write(*,*) 'Thank you.  Now please type a positive number which'
  215.       write(*,*) 'is close enough to zero (for use in identifying'
  216.       write(*,*) 'points on the boundary which should be on one of'
  217.       write(*,*) 'the contours).'
  218.       read(*,*) zero
  219.       write(*,*) 'Thank you.'
  220.       write(*,*) 'The parameter statements indicate that ',nlvl
  221.       write(*,*) 'contours will be drawn.  Please type the value of'
  222.       write(*,*) 'the dependent variable along the contour where that'
  223.       write(*,*) 'dependent variable is maximum.'
  224.       read(*,*) to
  225.       write(*,*) 'Thank you.  Please type the value of the dependent'
  226.       write(*,*) 'variable along the contour where that dependent'
  227.       write(*,*) 'variable is minimum.'
  228.       read(*,*) from
  229.       by = (to - from)/float(nlvl-1)
  230.       write(*,*) 'Thank you.  I compute the contour interval as',by
  231.       write(*,*) 'If something is wrong, please type 0 to stop the'
  232.       write(*,*) 'program.  Then edit the parameter statements.'
  233.       write(*,*) 'Otherwise, type 1 to continue with the program.'
  234.       read(*,*) iflag
  235.       if(iflag .eq. 0) stop
  236. c
  237. c     open files with variable file names to save wear and tear on user
  238.       open(9,file=incont)
  239.       open(10,file=outcon)
  240.       open(11,file=rdnbot)
  241.       open(12,file=outbnd)
  242.       open(13,file=labels)
  243. c
  244. c     read in the arrays x,z,h,nbot
  245.       read(9,*) (h(j),j=1,npts)
  246.       read(9,*) (x(j),j=1,npts)
  247.       read(9,*) (z(j),j=1,npts)
  248.       read(11,*) (nbot(jcol),jcol=0,lbot)
  249. c
  250. c     find the values of the dependent variable along the contours
  251.       write(*,*) 'values of the dependent variable along the contours'
  252.       do 100 i=1,nlvl
  253.            value(i) = from + (i-1)*by
  254.            len(i) = 0
  255.            write(*,*) value(i)
  256.   100 continue
  257.       return
  258.       end
  259. c
  260.       subroutine scan (h,x,z,nbot,value,len,ncell,xcont,zcont,npts,
  261.      $     nel,nlvl,lbot,mxlen,zero)
  262.       integer npts, nel,nlvl,lbot,mxlen
  263.       real h(npts), x(npts), z(npts), xcont(nlvl, mxlen),
  264.      $     zcont(nlvl,mxlen), value(nlvl)
  265.       integer nbot(0:lbot), len(nlvl), ncell(nlvl,nel)
  266. c     for each contour, scan all points in the array, looking for
  267. c     adjacent points where the dependent variable is bigger than
  268. c     the contour value at one point and smaller than the contour
  269. c     value at the other point.  Let the current point be a candidate
  270. c     for the upper lefthand corner of some element.  If the lower
  271. c     righthand corner and the upper righthand corner of that element
  272. c     exist, check for contour values between the upper lefthand corner
  273. c     and the lower lefthand corner, between the upper lefthand corner
  274. c     and the upper righthand corner, between the lower lefthand corner
  275. c     and the lower righthand corner, and between the upper righthand
  276. c     corner and the lower righthand corner.  If the lower righthand
  277. c     corner or the upper righthand corner does not exist, any needed
  278. c     contour points have already been found. Keep track of how many
  279. c     points are found for each element.  The appropriate element of
  280. c     ncell is the number of contour points in this cell for this
  281. c     contour level.  i counts the contour levels and icell counts the
  282. c     points in the contour level.
  283.       do 500 i=1,nlvl
  284.            val = value(i)
  285.            icell = 0
  286.            do 400 jcol = 0,lbot-2
  287.                 nbottm = nbot(jcol)
  288.                 nxtbot = nbot(jcol+1)
  289.                 lcol = nxtbot - nbottm
  290.                 ntop = nxtbot - 1
  291.                 ntopn = nbot(jcol+2) - 1
  292.                 do 300 j = nbottm + 1, ntop
  293. c                 do lower righthand corner and upper righthand corner
  294. c                 of this element exist?  If not, go to the next
  295. c                 element.
  296.                   jp = j+lcol
  297.                   if (jp .le. ntopn .and. jp-1 .le. ntopn) then
  298.                      icell = icell + 1
  299.                      ncell(i,icell) = 0
  300. c                    check for a contour point between the upper
  301. c                    lefthand corner and the lower lefthand corner
  302.                      if ( (val .lt. h(j) .and. val .gt. h(j-1)) .or.
  303.      $                    (val .gt. h(j) .and. val .lt. h(j-1)) ) then
  304. c                         the contour passes through a point between
  305. c                         points j and j-1.  Find it by linear
  306. c                         interpolation.
  307.                           frac = (h(j) - val) /(h(j) - h(j-1))
  308.                           dx = x(j) - x(j-1)
  309.                           dz = z(j) - z(j-1)
  310.                           xval = x(j) - frac*dx
  311.                           zval = z(j) - frac*dz
  312. c                         xval,zval is the point to plot.
  313.                           len(i) = len(i) + 1
  314.                           xcont(i,len(i)) = xval
  315.                           zcont(i,len(i)) = zval
  316.                           ncell(i,icell) = ncell(i,icell) + 1
  317.                      end if
  318.                      jp = j + lcol
  319. c                    check for a contour point between the upper
  320. c                    lefthand corner and the upper righthand corner.
  321.                      if ( (val .lt. h(j) .and. val .gt. h(jp)) .or.
  322.      $                      (val .gt. h(j) .and. val .lt. h(jp)) ) then
  323. c                           the contour passes through a point between
  324. c                           points j and jp.  Find it by linear
  325. c                           interpolation.
  326.                             frac = (h(j) - val) /(h(j) - h(jp))
  327.                             dx = x(j) - x(jp)
  328.                             dz = z(j) - z(jp)
  329.                             xval = x(j) - frac*dx
  330.                             zval = z(j) - frac*dz
  331. c                           xval,zval is the point to plot.
  332.                             len(i) = len(i) + 1
  333.                             xcont(i,len(i)) = xval
  334.                             zcont(i,len(i)) = zval
  335.                             ncell(i,icell) = ncell(i,icell) + 1
  336.                      end if
  337. c                    check for a contour point between the upper
  338. c                    righthand corner and the lower righthand corner.
  339.                      if ( (val .lt. h(jp) .and. val .gt. h(jp-1)) .or.
  340.      $                    (val .gt. h(jp) .and. val .lt. h(jp-1)) ) then
  341. c                         the contour passes through a point between
  342. c                         points jp and jp-1.  Find it by linear
  343. c                         interpolation.
  344.                           frac = (h(jp) - val) /(h(jp) - h(jp-1))
  345.                           dx = x(jp) - x(jp-1)
  346.                           dz = z(jp) - z(jp-1)
  347.                           xval = x(jp) - frac*dx
  348.                           zval = z(jp) - frac*dz
  349. c                         xval,zval is the point to plot.
  350.                           len(i) = len(i) + 1
  351.                           xcont(i,len(i)) = xval
  352.                           zcont(i,len(i)) = zval
  353.                           ncell(i,icell) = ncell(i,icell) + 1
  354.                      end if
  355. c                    check for a contour point between the lower
  356. c                    lefthand corner and the lower righthand corner.
  357.                      if ( (val .lt. h(j-1) .and. val .gt. h(jp-1)) .or.
  358.      $                      (val .gt. h(j-1) .and. val .lt. h(jp-1)) )
  359.      $                      then
  360. c                           the contour passes through a point between
  361. c                           points j-1 and jp-1.  Find it by linear
  362. c                           interpolation.
  363.                             frac = (h(j-1) - val) /(h(j-1) - h(jp-1))
  364.                             dx = x(j-1) - x(jp-1)
  365.                             dz = z(j-1) - z(jp-1)
  366.                             xval = x(j-1) - frac*dx
  367.                             zval = z(j-1) - frac*dz
  368. c                           xval,zval is the point to plot.
  369.                             len(i) = len(i) + 1
  370.                             xcont(i,len(i)) = xval
  371.                             zcont(i,len(i)) = zval
  372.                             ncell(i,icell) = ncell(i,icell) + 1
  373.                      end if
  374.                      if (abs(val - h(j)) .le. zero) then
  375. c                         check for a contour point at j
  376. c                         (upper lefthand corner)
  377. c                         check for val being equal to h(j)
  378.                           len(i) = len(i) + 1
  379.                           xcont(i,len(i)) = x(j)
  380.                           zcont(i,len(i)) = z(j)
  381.                           ncell(i,icell) = ncell(i,icell) + 1
  382.                      end if
  383.                      if (abs(val - h(j-1)) .le. zero) then
  384. c                         check for a contour point at j-1
  385. c                         (lower lefthand corner)
  386. c                         check for val being equal to h(j-1)
  387.                           len(i) = len(i) + 1
  388.                           xcont(i,len(i)) = x(j-1)
  389.                           zcont(i,len(i)) = z(j-1)
  390.                           ncell(i,icell) = ncell(i,icell) + 1
  391.                      end if
  392.                      if (abs(val - h(jp)) .le. zero) then
  393. c                         check for a contour point at jp
  394. c                         (upper righthand corner)
  395. c                         check for val being equal to h(jp)
  396.                           len(i) = len(i) + 1
  397.                           xcont(i,len(i)) = x(jp)
  398.                           zcont(i,len(i)) = z(jp)
  399.                           ncell(i,icell) = ncell(i,icell) + 1
  400.                      end if
  401.                      if (abs(val - h(jp-1)) .le. zero) then
  402. c                         check for a contour point at jp-1
  403. c                         (lower righthand corner)
  404. c                         check for val being equal to h(jp-1)
  405.                           len(i) = len(i) + 1
  406.                           xcont(i,len(i)) = x(jp-1)
  407.                           zcont(i,len(i)) = z(jp-1)
  408.                           ncell(i,icell) = ncell(i,icell) + 1
  409.                      end if
  410.                   end if
  411.   300           continue
  412.   400      continue
  413.   500 continue
  414.       return
  415.       end
  416. c
  417.       subroutine wrcon(ncell,xcont,zcont,last,nlvl,nel,mxlen)
  418.       integer nlvl,nel,mxlen
  419.       real xcont(nlvl,mxlen), zcont(nlvl,mxlen)
  420.       integer ncell(nlvl,nel), last(nlvl)
  421. c     write information needed by gnuplot to file in format appropriate
  422. c     for gnuplot.  Each blank line causes gnuplot to lift the pen.
  423. c     Write the contour points for each cell with the order and
  424. c     repetition required for Sundberg's algorithm.  Then lift the
  425. c     pen and go to the next cell.
  426.       do 700 i=1,nlvl
  427.            last(i) = 0
  428.            do 650 icell = 1,nel
  429.            if (ncell(i,icell) .eq. 2) then
  430.                 write(10,*) xcont(i,1+last(i)),
  431.      $                      zcont(i,1+last(i))
  432.                 write(10,*) xcont(i,2+last(i)),
  433.      $                      zcont(i,2+last(i))
  434.                 write(10,*)
  435.                 last(i) = last(i) + 2
  436.            else if (ncell(i,icell) .eq. 4) then
  437.                 write(10,*) xcont(i,1+last(i)),
  438.      $                      zcont(i,1+last(i))
  439.                 write(10,*) xcont(i,2+last(i)),
  440.      $                      zcont(i,2+last(i))
  441.                 write(10,*) xcont(i,3+last(i)),
  442.      $                      zcont(i,3+last(i))
  443.                 write(10,*) xcont(i,4+last(i)),
  444.      $                      zcont(i,4+last(i))
  445.                 write(10,*) xcont(i,1+last(i)),
  446.      $                      zcont(i,1+last(i))
  447.                 write(10,*) xcont(i,3+last(i)),
  448.      $                      zcont(i,3+last(i))
  449.                 write(10,*) xcont(i,2+last(i)),
  450.      $                      zcont(i,2+last(i))
  451.                 write(10,*) xcont(i,4+last(i)),
  452.      $                      zcont(i,4+last(i))
  453.                 write(10,*)
  454.                 last(i) = last(i) + 4
  455.            end if
  456.   650      continue
  457.   700 continue
  458.       return
  459.       end
  460. c
  461.       subroutine wrbdy(x,z,nbot,npts,lbot)
  462.       integer npts, lbot
  463.       real x(npts), z(npts)
  464.       integer nbot(0:lbot)
  465. c     Write the coordinates of the points on the boundary to another
  466. c     file.
  467.       do 900 j=1,nbot(1)-1
  468.            write(12,*) x(j),z(j)
  469.   900 continue
  470.       do 1000 jj=1,lbot
  471.            j=nbot(jj)-1
  472.            write(12,*) x(j),z(j)
  473.  1000 continue
  474.       do 1100 j=nbot(lbot)-1,nbot(lbot-1),-1
  475.            write(12,*) x(j),z(j)
  476.  1100 continue
  477.       do 1200 jj=lbot-1,0,-1
  478.            j=nbot(jj)
  479.            write(12,*) x(j),z(j)
  480.  1200 continue
  481.       return
  482.       end
  483. c
  484.       subroutine wrlabl(xcont,zcont,value,len,nlvl,mxlen)
  485.       integer nlvl, mxlen
  486.       real xcont(nlvl,mxlen), zcont(nlvl,mxlen), value(nlvl)
  487.       integer len(nlvl)
  488. c     Write the labels for the contours to a file for use by gnuplot
  489.       write(13,*) 'set parametric'
  490.       write(13,*) 'set data style lines'
  491.       ii = 10
  492. c     (allow labels 1 through 9 for other information)
  493. c     If necessary, remove some of the labels manually or add new ones
  494. c     Modify the labels manually to remove blank spaces and superfluous
  495. c     zeros.
  496.       do 1400 i=1,nlvl
  497.            if (len(i) .ne. 0) then
  498. c               write the label once, in the middle of the contour
  499.                 write(13,1300) ii,value(i), xcont(i,len(i)/2 + 1),
  500.      $                         zcont(i,len(i)/2 + 1)
  501.  1300           format('set label ',i3,' "',f10.3,'" at ',f10.3,
  502.      $                 ',',f10.3,' center')
  503.                 ii = ii + 1
  504.            end if
  505.  1400 continue
  506.       return
  507.       end
  508.